Project cuối kì
2024-7-7
EDA + Preprocessing
data <- read.csv('data/SeoulBikeData.csv',fileEncoding = "latin1",check.names = F)
glimpse(data)## Rows: 8,760
## Columns: 14
## $ Date <chr> "01/12/2017", "01/12/2017", "01/12/2017", …
## $ `Rented Bike Count` <int> 254, 204, 173, 107, 78, 100, 181, 460, 930…
## $ Hour <int> 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, …
## $ `Temperature(°C)` <dbl> -5.2, -5.5, -6.0, -6.2, -6.0, -6.4, -6.6, …
## $ `Humidity(%)` <int> 37, 38, 39, 40, 36, 37, 35, 38, 37, 27, 24…
## $ `Wind speed (m/s)` <dbl> 2.2, 0.8, 1.0, 0.9, 2.3, 1.5, 1.3, 0.9, 1.…
## $ `Visibility (10m)` <int> 2000, 2000, 2000, 2000, 2000, 2000, 2000, …
## $ `Dew point temperature(°C)` <dbl> -17.6, -17.6, -17.7, -17.6, -18.6, -18.7, …
## $ `Solar Radiation (MJ/m2)` <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, …
## $ `Rainfall(mm)` <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ `Snowfall (cm)` <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ Seasons <chr> "Winter", "Winter", "Winter", "Winter", "W…
## $ Holiday <chr> "No Holiday", "No Holiday", "No Holiday", …
## $ `Functioning Day` <chr> "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", …data <- data |> janitor::clean_names()
data <- data %>%
filter(functioning_day == "Yes") |> select(-functioning_day)
data <- data %>%
drop_na()
data <- data %>%
mutate(holiday = ifelse(holiday == "No Holiday", 0, 1))
data <- data %>%
mutate(time = case_when(
hour >=4 & hour <11 ~ "sang",
hour >=11 & hour <13 ~ "trua",
hour >=13 & hour <18 ~ "chieu",
.default = "toi"
))
data$seasons <- factor(data$seasons)
data$holiday <- factor(data$holiday)
data$hour <- factor(data$hour)Visualization
data |> select(-c(date,hour,seasons,holiday,time)) |> gather(variable, value) |> group_by(variable)|> summarise(mean = round(mean(value),4),sd = sd(value),median =median(value))## # A tibble: 9 × 4
## variable mean sd median
## <chr> <dbl> <dbl> <dbl>
## 1 dew_point_temperature_c 3.94 13.2 4.7
## 2 humidity_percent 58.1 20.5 57
## 3 rainfall_mm 0.149 1.13 0
## 4 rented_bike_count 729. 642. 542
## 5 snowfall_cm 0.0777 0.444 0
## 6 solar_radiation_mj_m2 0.568 0.868 0.01
## 7 temperature_c 12.8 12.1 13.5
## 8 visibility_10m 1434. 609. 1690
## 9 wind_speed_m_s 1.73 1.03 1.5data_m <-data |> select(-c(date)) |> gather(variable, value,-c(rented_bike_count,hour,seasons,holiday,time))
ggplot(data_m,aes(x = value, y = after_stat(density))) +
facet_wrap(~variable, scales = "free") +
geom_histogram(aes(y = after_stat(density)), color = "white", fill="lightblue", alpha=1, position="identity", add_density = TRUE) + geom_density() + theme_minimal()
Tổng quan về dữ liệu:
Dew point temperature: Dữ liệu phân bố từ -30°C tới 25°C. Trong đó, phần lớn giá trị nằm trong khoảng nhiệt độ từ 0°C đến 25°C.
Humidity percent: Biểu đồ histogram cho thấy phân bố độ ẩm phần trăm từ 0% đến 100%. Các giá trị độ ẩm thường xuất hiện nhất rơi vào khoảng từ 40% đến 80%.
Rainfall: Lượng mưa chủ yếu tập trung ở gần giá trị 0. Dựa vào biểu đồ, có thể thấy các giá trị mưa lớn là rất hiếm trong tập dữ liệu.
Snowfall: Dữ liệu phân bố từ 0cm tới 9cm. Tuy nhiên, phần lớn dữ liệu tập trung ở giá trị 0cm.
Solar radiation: Phân bố bức xạ mặt trời tập trung chủ yếu ở giá trị rất thấp gần 0 MJ/m².
Temperature: Dữ liệu nhiệt độ phân bố từ -17°C tới 40°C. Trong đó, nhiệt độ phân bố chủ yếu từ 0°C tới 25°C. Giá trị 20°C có mật độ cao nhất.
Visibility: Biểu đồ cho thấy phần lớn các giá trị tầm nhìn tập trung ở khoảng 2000 đơn vị. Ngoài một số ít trường hợp, dữ liệu cho thấy tần suất xuất hiện của các giá trị tầm nhìn thấp hơn (dưới 2000 đơn vị) là rất ít.
Wind speed: Biểu đồ cho thấy tốc độ gió thường tập trung nhiều ở khoảng từ 0 đến 2m/s, với đỉnh tại khoảng 1m/s.
ggplot(data_m,aes(x = value,color = seasons)) +
facet_wrap(~variable, scales = "free") +
geom_boxplot() + theme_minimal()
Nhận xét về các boxplot:
Dew point temperature: Các điểm ngưng sương cao nhất ở mùa hè, thấp nhất ở mùa đông, các mùa còn lại có trung vị là trung bình của 2 mùa hè và đông
Humidity percent: Thấp nhất tại mùa đông và độ ẩm cao tương đương nhau ở các mùa còn lại
Rainfall: Lưu lượng mưa có giá trị có thể quan sát cao nhất tại mùa hè, xuân, thu còn mùa đông lượng mưa không có quan sát
Snowfall: Độ dày của tuyết chỉ có thể quan sát tại mùa đông
Solar radiation: Bức xạ mặt trời cao tại mùa hè và thấp ở mùa đông
Temperature: Nhiệt độ cao tại mùa hè và thấp tại mùa đông, 2 mùa còn lại tương đương nhau
Visbility: Cao nhất tại mùa thu và thấp nhất tại mùa xuân, mùa đông và hè tương đương
Wind speed: Tốc độ gió gần như không có khác biệt đáng kể giữa các mùa.
data_by_day <- data %>% mutate(date = as.Date(date, "%d/%m/%Y")) %>% group_by(date) %>% summarise(rented_bike_count = sum(rented_bike_count), seasons = first(seasons))
ggplot(data = data_by_day, aes(x=date, y=rented_bike_count, colour = seasons)) + geom_line(alpha=.3) + geom_smooth(span = 0.2, se=F) + theme_minimal()
Nhận xét:
Vào mùa đông, số lượng thuê xe đạp là thấp nhất.
Số lượng thuê xe đạp cao nhất ở thời điểm cuối mùa xuân và đầu mùa hè.
Số lượng thuê xe đạp có sự tăng, giảm không ổn định. Điều này có thể do tác động của những yếu tố như ngày nghỉ, lượng mưa, độ ẩm,….
ggplot(data_m,aes(y = rented_bike_count,x = value, color = time)) +
facet_wrap(~variable, scales = "free") +
geom_point( alpha = 0.3, shape = 16)+
stat_smooth(se=F) 
Nhận xét: Nhìn chung, mối quan hệ giữa các biến độc lập và biến phụ thuộc tương đối ổn định qua các thời điểm trong ngày (trừ biến Solar radition).
ggplot(data_m,aes(y = rented_bike_count,x = value, color = seasons)) +
facet_wrap(~variable, scales = "free") +
geom_point( alpha = 0.3, shape = 16)+
stat_smooth(se=F)
Nhận xét: Có thể thấy, mối quan hệ giữa các biến độc lập và biến phụ thuộc có sự thay đổi khá nhiều qua từng mùa.
** Vì vậy, chúng ta cần xây dựng mô hình cho từng mùa để tránh sự biến động của dữ liệu, tăng sự ổn định của mô hình hồi quy và giảm độ phức tạp trong phân tích. **
** Ngoài ra, dựa vào các biểu đồ, có thể thấy mối quan hệ giữa các biến độc lập và biến phụ tuân theo hàm mũ, do đó, chúng ta nên áp dụng phương pháp hồi quy poisson. **
ggplot(data,aes(y = humidity_percent,x = hour, color = rented_bike_count)) +
geom_jitter(alpha = 0.5, shape = 16) 
Nhận xét:
Từ 0h đến 9h và từ 19h đến 23h: Có thể thấy trong hai khung giờ này, khi độ ẩm nằm trong khoảng từ 50% đến 80%, số lượng khách hàng thuê xe đạp nhiều hơn. Điều này có thể là do trong hai khung giờ đó, nhiệt độ thường khá mát mẻ hoặc lạnh, và độ ẩm từ 50% đến 80% là điều kiện lý tưởng để đi xe đạp.
Từ 9h đến 19h: Ở khung giờ này, thời tiết đang nóng dần lên nên khách hàng có nhu cầu thuê xe đạp nhiều hơn khi độ ẩm nằm trong khoảng từ 25% tới 60%.
ggplot(data,aes(y = temperature_c,x = hour, color = rented_bike_count)) +
geom_jitter(alpha = 0.5, shape = 16) 
Nhận xét:
Số lượng xe đạp được thuê có xu hướng tăng vào các giờ sáng sớm (khoảng 7-9 giờ) và buổi chiều tối (khoảng 17-19 giờ).
Khi nhiệt độ tăng, đặc biệt là từ khoảng 10°C đến 30°C, số lượng xe đạp được thuê cũng tăng lên.
Ở các nhiệt độ cực đoan (dưới 0°C và trên 30°C), số lượng xe đạp được thuê giảm đi, điều này có thể do điều kiện thời tiết không thuận lợi cho việc đi xe đạp.
ggplot(data,aes(x = rented_bike_count,y = time)) +
stat_halfeye( position = "dodge" )+
facet_wrap(~seasons)
Nhận xét:
Mùa thu: Không có sự khác biệt rõ rệt giữa các thời điểm trong ngày. Số lượng xe đạp được thuê thấp và phân bố tương đối đồng đều.
Mùa xuân: Buổi sáng và buổi tối có số lượng xe đạp được thuê thấp, với phân bố khá nhỏ. Buổi chiều có số lượng xe đạp được thuê cao nhất, điều này có thể do người dân tham gia nhiều hoạt động ngoài trời hơn sau giờ làm việc và khi thời tiết ấm dần lên trong ngày. Buổi trưa có số lượng xe đạp được thuê cao hơn so với sáng và tối, nhưng vẫn thấp hơn buổi chiều.
Mùa hè: Buổi trưa có số lượng xe đạp được thuê cao nhất. Buổi tối và buổi chiều có số lượng xe đạp được thuê thấp hơn, có thể do nhiệt độ buổi chiều và buổi tối không thuận lợi hoặc người dân ít có nhu cầu thuê xe vào các thời điểm này.
Mùa đông: Số lượng xe đạp được thuê rất ít vào tất cả các thời điểm trong ngày, có thể do điều kiện thời tiết khắc nghiệt, lạnh giá làm giảm nhu cầu sử dụng xe đạp.
Nhận xét: Số lượng thuê xe đạp thay đổi theo từng khung thời gian.
Nhận xét:
Vào mùa Thu và Xuân, phân phối số lượng xe đạp thuê có mô hình tương tự với mật độ cao hơn ở các giá trị thấp và đuôi dài mở rộng đến các giá trị cao hơn.
Mùa Hạ có phân phối tương đối hẹp hơn, cho thấy sự biến đổi ít hơn trong việc thuê xe đạp, với số lượng thuê vào các ngày nghỉ lễ cao hơn một chút so với các ngày không phải nghỉ lễ.
Mùa Đông cho thấy số lượng xe đạp thuê thấp nhất với sự giảm rõ rệt so với các mùa khác. Phân phối cho các ngày nghỉ lễ và không nghỉ lễ khá giống nhau vào mùa Đông.
Preprocessing
data.clean1 <- data |> group_by(seasons,hour) |>
mutate(l = quantile(rented_bike_count,0.25) - 1.5*IQR(rented_bike_count),
h = quantile(rented_bike_count,0.75) + 1.5*IQR(rented_bike_count)) |>
filter(rented_bike_count >= l & rented_bike_count <= h) |> select(-c(h,l)) |> ungroup()dis_mahalanobis <- function(X){
return(mahalanobis(X, colMeans(X),diag(1e-9, ncol(X))+ cov(X)))
}
index <- data.clean1 |> select(-rented_bike_count)|> mutate(id= row_number()) |>
group_by(seasons,hour) |>
mutate(mahalanobis_dist = dis_mahalanobis(across(where(is.numeric))) ) |>
filter( mahalanobis_dist <= qchisq(0.95, ncol(across(where(is.numeric))))) |>
pull(id)
data.clean <- data.clean1[index,] |> select(-c(date))
nrow(data.clean)## [1] 7872Testing
Nhận định ban đầu
Khi ta xét đến các bài toán về việc quản lý các dịch vụ công, ta cần quan tâm tới những ngày mà số lượng người sử dụng dịch vụ tăng mạnh. Thông thường ta có thể thấy vào những ngày nghỉ, ngày lễ số lượng người ra đường vui chơi sẽ nhiều hơn so với ngày thường. Từ giả định đó ta sẽ kiểm tra xem việc ngày đó có phải ngày nghỉ hay không liệu có ảnh hưởng đến số lượng xe đạp được thuê hay không.
AB test cho holiday
Ta xem xét câu hỏi: “Liệu số lượng xe đạp được thuê vào các ngày nghỉ, ngày lễ có tăng đột biến hay không để từ đó thành phố có thể đề ra các kế hoạch ứng biến vào các ngày đặc biệt này”.
Trước hết ta có thể thấy rõ rằng phần trăm ngày nghỉ so với các ngày trong năm là rất thấp (ví dụ trong 1 tháng 4 tuần chỉ có 4 ngày Chủ Nhật so với 24 ngày bình thường), tuy nhiên không thể bỏ qua các trường hợp người ta sẽ gia tăng đột ngột việc thuê xe để giải trí, vui chơi vào các ngày nghỉ và ngày lễ.
Do đó ta cần kiểm định cho hai loại ngày này.
data <- data.clean
data$holiday <- as.factor(data$holiday)
data |> group_by(holiday) |> summarise(n = n(), mean = mean(rented_bike_count), sd = sd(rented_bike_count))## # A tibble: 2 × 4
## holiday n mean sd
## <fct> <int> <dbl> <dbl>
## 1 0 7485 763. 643.
## 2 1 387 533. 562.Kết quả cho thấy có 7485 lượt thuê xe vào ngày thường và 387 lượt thuê xe vào ngày nghỉ. Trung bình lượng thuê xe với ngày thường là cao hơn ngày nghỉ. Độ lệch chuẩn của ngày thường cũng lớn hơn, cho thấy độ biến động lượng thuê xe giữa các ngày thường là lớn hơn so với ngày nghỉ
Biểu đồ violin dưới đây giúp ta khẳng định các nhận định trên, đồng thời cung cấp thêm thông tin về phân phối lượng thuê xe đạp, ở đây, dữ liệu của cả hai ngày đều cho thấy phân phối bất đối xứng và lệch phải của lượng thuê xe.
ggplot(data, aes(x = holiday, y = rented_bike_count, fill = holiday)) +
geom_violin() +
geom_boxplot(width = 0.15) +
labs(x = "Holiday", y = "Rented bike count") +
theme_bw() +
theme(legend.position = "none")
Thông qua bảng tổng hợp và biểu đổ violin, một giả định có thể được đưa ra: “số lượng thuê xe của ngày thường là nhiều hơn ngày nghỉ”. Do đó, ta cần kiểm chứng giả thuyết và đối thuyết sau:
\(H_0: \mu_1 = \mu_0\)
\(H_1: \mu_1 < \mu_0\)
Nếu \(H_0\) đúng tức sự khác biệt giữa số lượng xe đạp được thuê giữa ngày nghỉ và ngày thường chỉ là kết quả của sự ngẫu nhiên (không có ý nghĩa thống kê). Để kiểm định giả thuyết này ta sử dụng kiểm định hoán vị, và p-value sẽ được tính cho kiểm định bên trái.
perm_fun <- function(x, n1, n0, R) {
n <- n1 + n0
mean_diff <- numeric(R)
for (i in 1:R){
idx_1 <- sample(x = 1:n, size = n1)
idx_0 <- setdiff(x = 1:n, y = idx_1)
mean_diff[i] <- mean(x[idx_1]) - mean(x[idx_0])
}
return(mean_diff)
}
set.seed(42)
diff_mean_perm <- perm_fun(data$rented_bike_count, n1 = 387, n0 = 7485, R = 10000)
ggplot(data = tibble(perm_diffs = diff_mean_perm), aes(x = perm_diffs)) +
geom_histogram(bins = 10, fill = "gray80", color = "black") +
labs(x = "Rented bike count differences", y = "Frequency") +
theme_bw()
Giá trị p-value
mean_0 <- mean(data$rented_bike_count[data$holiday == '0'])
mean_1 <- mean(data$rented_bike_count[data$holiday == '1'])
mean(diff_mean_perm < (mean_1 - mean_0))## [1] 0Với mức ý nghĩa \(\alpha = 0.05\), kết quả cho thấy sự khác biệt giữa số lượng xe đạp được thuê trong ngày nghỉ và ngày thường không phải do sự ngẫu nhiên (có ý nghĩa thống kê).
Từ kết luận trên ta có thể thấy rằng, số lượng xe đạp được thuê vào ngày nghỉ sẽ không ảnh hưởng nhiều (không tăng đột biến) và vẫn ít hơn ngày thường. Tuy nhiên đây không phải yếu tố ảnh hưởng duy nhất đến việc thuê xe đạp của người dân
ANOVA test cho seasons
Tiếp đến, cùng câu hỏi như trên nhưng lần này ta lại xét đến các mùa trong năm: “số lượng thuê xe đạp có bị ảnh hưởng bởi các mùa hay không?”
Việc các mùa trong năm có thể ảnh hưởng đến số lượng thuê xe đạp rất nhiều. Ta có thể phỏng đoán rằng vào mùa đông thời tiết giá lạnh, người ta sẽ lựa chọn các loại phương tiện khác để giữ ấm thay vì chọn đạp xe dưới thời tiết khắc nghiệt kể trên. Hoặc vào mùa hè một lượng lớn du khách nước ngoài đến du lịch và trải nghiệm đạp xe ở Seoul cũng góp phần làm thay đổi số lượng thuê xe so với các mùa khác.
Do đó ta sẽ kiểm định cho cả bốn mùa trong năm.
data |> group_by(seasons) |> summarise(n = n(), mean = mean(rented_bike_count), sd = sd(rented_bike_count))## # A tibble: 4 × 4
## seasons n mean sd
## <fct> <int> <dbl> <dbl>
## 1 Autumn 1756 971. 604.
## 2 Spring 2057 755. 612.
## 3 Summer 2023 1085. 673.
## 4 Winter 2036 228. 148.Đúng như chúng ta dự đoán trung bình số lượng xe đạp được thuê vào mùa đông là thấp nhất và vào mùa hè là cao nhất. Ta có thể thấy trung bình số lượng xe đạp được thuê của 4 mùa là khác nhau (về mặt giá trị).
Trong khi đó, độ lệch chuẩn là cho thấy độ biến động trong số lượng xe đạp được thuê là khác biệt giữa các nhóm.
Biểu đồ violin làm rõ hơn cho nhận xét ở trên.
ggplot(data, aes(x = seasons, y = rented_bike_count, fill = seasons)) +
geom_violin() +
geom_boxplot(width = 0.15) +
labs(x = "Seasons", y = "Rented bike count") +
theme_bw() +
theme(legend.position = "none")
Ta đặt giả thuyết như sau:
\(H_0: \mu_{\text{autumn}} = \mu_{\text{spring}} = \mu_{\text{summer}} = \mu_{\text{winter}}\)
\(H_1:\) Ít nhất có một trung bình là khác với những cái còn lại.
Để kiểm định cho giả thuyết trên ta thực hiện Permutation ANOVA. Permutation ANOVA là một phương pháp phân tích ANOVA khác không phụ thuộc vào một số giả định như là: dữ liệu trong mỗi nhóm phải tuân theo phân phối chuẩn, giả định về sự đồng nhất phương sai,…
set.seed(42)
out_aov_1 <- aovp(formula = rented_bike_count ~ seasons, data = data, perm = "Prob")## [1] "Settings: unique SS "summary(out_aov_1)## Component 1 :
## Df R Sum Sq R Mean Sq Iter Pr(Prob)
## seasons 3 867441627 289147209 5000 < 2.2e-16 ***
## Residuals 7868 2369546013 301162
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Giá trị p-value được cho cung cấp bởi cột Pr(Prob) là 2.2e-16. Số lần lấy mẫu lặp lại là 5000, được cung cấp bởi cột Iter.
Với mức ý nghĩa \(\alpha = 0.05\), kết quả cho thấy sự khác biệt giữa số lượng xe đạp được thuê vào từng mùa là có ý nghĩa thống kê (có ít nhất một trung bình là khác với các cái còn lại).
Từ kết luận trên, ta có thể thấy các mùa trong năm có ảnh hưởng đến số lượng thuê xe đạp. Thành phố có thể đưa ra các kế hoạch như là tăng cường xe đạp vào mùa hè để đáp ứng đủ nhu cầu và giảm bớt số lượng xe đạp vào mùa đông (có thể đưa về kho để tránh tình trạng hư hại do thời tiết).
ANOVA test kết hợp
Từ các dữ kiện ở trên ta có thể thực hiện ANOVA cho cả mùa và ngày nghỉ. Ta sẽ xem xét thêm rằng các ngày nghỉ trong một mùa có ảnh hưởng tới số lượng thuê xe đạp hay không? Ví dụ vào mùa hè các học sinh được nghỉ học tuy nhiên bố mẹ chúng vẫn đi làm, chỉ cho đến cuối tuần mới có thời gian để đưa con của họ đi chơi, số lượng thuê xe đạp có thể sẽ tăng vào những ngày như vậy.
data |> group_by(seasons, holiday) |> summarise(n = n(), mean = mean(rented_bike_count), sd = sd(rented_bike_count))## # A tibble: 8 × 5
## # Groups: seasons [4]
## seasons holiday n mean sd
## <fct> <fct> <int> <dbl> <dbl>
## 1 Autumn 0 1663 975. 605.
## 2 Autumn 1 93 914. 581.
## 3 Spring 0 1995 758. 612.
## 4 Spring 1 62 680. 605.
## 5 Summer 0 1975 1087. 675.
## 6 Summer 1 48 1022. 564.
## 7 Winter 0 1852 235. 150.
## 8 Winter 1 184 163. 106.Có vẻ như phỏng đoán của chúng ta lại chính xác khi mà số ngày nghỉ, ngày lễ trong mùa hè rất thấp so với ngày thường tuy nhiên trung bình lượng thuê xe đạp lại sấp sỉ với ngày thường. Điều này cũng đúng với trường hợp mùa thu (có thể do mùa này không khí dễ chịu, cảnh sắc tươi đẹp nên người ta thường dành cuối tuần để đạp xe thư giãn)
ggplot(data, aes(x = seasons, y = rented_bike_count, fill = holiday)) +
geom_violin() +
labs(x = "Seasons", y = "Rented bike count", fill = "Holiday") +
theme_bw() +
theme(legend.position = "bottom") +
facet_wrap(~seasons, scales = "free") +
geom_boxplot(width = 0.15, position = position_dodge(width = 0.89)) 
Để kiểm tra xem việc này là do ngẫu nhiên hay không, ta sẽ thực hiện Permutation ANOVA.
set.seed(42)
out_aov_2 <- aovp(formula = rented_bike_count ~ seasons*holiday, data = data, perm = "Prob")## [1] "Settings: unique SS "summary(out_aov_2)## Component 1 :
## Df R Sum Sq R Mean Sq Iter Pr(Prob)
## seasons 3 190825671 63608557 5000 <2e-16 ***
## holiday 1 1371580 1371580 5000 0.0126 *
## seasons:holiday 3 12797 4266 51 1.0000
## Residuals 7864 2367786918 301092
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Ta có thể thấy rằng, như kết luận ở phần trước, khi đứng độc lập thì các mùa trong năm và loại ngày sẽ ảnh hưởng đến số lượng thuê xe đạp (với mức ý nghĩa 5%).
Ở kiểm định này ta cũng có thể thấy thêm được sự tương tác giữa ngày nghỉ và các mùa trong năm. Với mức ý nghĩa \(\alpha = 0.05\), kết quả (p-value = 1) cho thấy sự khác biệt giữa số lượng xe đạp được thuê vào các loại ngày theo từng mùa là không có ý nghĩa thống kê.
Do đó các ngày nghỉ, ngày lễ vẫn không phải là tác nhân chính ảnh hưởng đến việc thuê xe của người dân Hàn Quốc.
Mở rộng
Việc ngày nghỉ, ngày lễ chỉ chiếm một phần nhỏ trong năm do đó khó mà gây ảnh hưởng đến số lượng xe đạp của người dân. Ta tiếp tục xét đến yếu tố thời gian trong ngày, ta sẽ kiểm định xem thời gian trong ngày có ảnh hưởng đến số lượng thuê xe hay không
data$holiday <- as.factor(data$time)
data |> group_by(time) |> summarise(n = n(), mean = mean(rented_bike_count), sd = sd(rented_bike_count))## # A tibble: 4 × 4
## time n mean sd
## <chr> <int> <dbl> <dbl>
## 1 chieu 1626 927. 565.
## 2 sang 2311 510. 503.
## 3 toi 3282 845. 746.
## 4 trua 653 707. 375.Ta có thể thấy trung bình số lượng xe đạp được thuê vào buổi chiều (từ 13h đến 18h) là cao nhất
ggplot(data, aes(x = time, y = rented_bike_count, fill = time)) +
geom_violin() +
geom_boxplot(width = 0.15) +
labs(x = "Seasons", y = "Rented bike count") +
theme_bw() +
theme(legend.position = "none")
set.seed(42)
out_aov_3 <- aovp(formula = rented_bike_count ~ time, data = data, perm = "Prob")## [1] "Settings: unique SS "summary(out_aov_3)## Component 1 :
## Df R Sum Sq R Mean Sq Iter Pr(Prob)
## time1 3 214695712 71565237 5000 < 2.2e-16 ***
## Residuals 7868 3022291929 384125
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Với mức ý nghĩa \(\alpha = 0.05\), kết quả cho thấy sự khác biệt giữa số lượng xe đạp được thuê vào từng khung giờ là có ý nghĩa thống kê (có ít nhất một trung bình là khác với các cái còn lại).
Từ kết luận trên, ta có thể thấy các khung giờ trong ngày có ảnh hưởng đến số lượng thuê xe đạp. Thành phố có thể đưa ra các kế hoạch như là tăng cường xe đạp vào từng khung giờ cụ thể để đáp ứng đủ nhu cầu của người dân.
Tương tự như trên ta cũng kiểm tra xem các khung giờ theo từng mùa trong năm
data |> group_by(seasons, holiday) |> summarise(n = n(), mean = mean(rented_bike_count), sd = sd(rented_bike_count))## # A tibble: 16 × 5
## # Groups: seasons [4]
## seasons holiday n mean sd
## <fct> <fct> <int> <dbl> <dbl>
## 1 Autumn chieu 356 1276. 421.
## 2 Autumn sang 512 661. 548.
## 3 Autumn toi 750 1038. 665.
## 4 Autumn trua 138 974. 212.
## 5 Spring chieu 429 1044. 541.
## 6 Spring sang 594 491. 478.
## 7 Spring toi 862 789. 689.
## 8 Spring trua 172 777. 357.
## 9 Summer chieu 428 1120. 488.
## 10 Summer sang 599 730. 536.
## 11 Summer toi 821 1374. 762.
## 12 Summer trua 175 865. 275.
## 13 Winter chieu 413 306. 105.
## 14 Winter sang 606 184. 171.
## 15 Winter toi 849 218. 144.
## 16 Winter trua 168 253. 78.7Tất cả các mùa đều có trung bình số lượng xe đạp được thuê vào buổi chiều là cao nhất.
ggplot(data, aes(x = seasons, y = rented_bike_count, fill = time)) +
geom_violin() +
labs(x = "Seasons", y = "Rented bike count", fill = "Time") +
theme_bw() +
theme(legend.position = "bottom") +
facet_wrap(~seasons, scales = "free") +
geom_boxplot(width = 0.15, position = position_dodge(width = 0.89)) 
set.seed(42)
out_aov_4 <- aovp(formula = rented_bike_count ~ seasons*time, data = data, perm = "Prob")## [1] "Settings: unique SS "summary(out_aov_4)## Component 1 :
## Df R Sum Sq R Mean Sq Iter Pr(Prob)
## seasons 3 548011049 182670350 5000 < 2.2e-16 ***
## time1 3 217347833 72449278 5000 < 2.2e-16 ***
## seasons:time1 9 109074873 12119430 5000 < 2.2e-16 ***
## Residuals 7856 2048677419 260779
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Với mức ý nghĩa \(\alpha = 0.05\), kết quả cho thấy sự khác biệt giữa số lượng xe đạp được thuê vào từng khung giờ theo từng mùa là có ý nghĩa thống kê, cho thấy rằng tác động của mùa lên số lượng xe đạp thuê thay đổi tùy thuộc vào từng khung giờ.
Regression
Data Model
Xây dựng mô hình hồi quy tuyến tính với biến phụ thuộc là rented_bike_count, và biến giải thích là các biến còn lại. 2 hướng nhắm tới là Linear Models và Generalized Linear Models. Chia tập dữ liệu thành 2 phần để đánh giá.
data.model <- data.clean |> select(-time)
slit_test_train <- function(df,pro=0.8,seed =30){
set.seed(seed)
N <- nrow(df)
ind_train <- sample(1:N,size=floor(pro*N))
ind_test <- setdiff(1:N,ind_train)
data_train <- tibble( df[ind_train,])
data_test <- tibble( df[ind_test,])
return(list(data_train,data_test))
}
dummies <- dummyVars("~ .", data = data.model |> select(-rented_bike_count))
data.numeric <- data.frame(predict(dummies, newdata = data.model |> select(-rented_bike_count)))
data.numeric$rented_bike_count <- data.model$rented_bike_count
rs <- slit_test_train(data.numeric,0.8)
trainset_numeric <- as.data.frame(rs[1])
testset_numeric <- as.data.frame(rs[2])
rs <- slit_test_train(data.model,0.8)
trainset <- as.data.frame(rs[1])
testset <- as.data.frame(rs[2])Linear Regression
md.ln <- lm(rented_bike_count ~ .,
data = trainset)
summary(md.ln)##
## Call:
## lm(formula = rented_bike_count ~ ., data = trainset)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1352.97 -215.39 -11.04 200.77 1725.97
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.306e+03 1.230e+02 10.618 < 2e-16 ***
## hour1 -1.231e+02 3.211e+01 -3.835 0.000127 ***
## hour2 -2.399e+02 3.199e+01 -7.499 7.31e-14 ***
## hour3 -3.060e+02 3.191e+01 -9.589 < 2e-16 ***
## hour4 -3.816e+02 3.182e+01 -11.993 < 2e-16 ***
## hour5 -3.580e+02 3.229e+01 -11.084 < 2e-16 ***
## hour6 -1.970e+02 3.217e+01 -6.126 9.58e-10 ***
## hour7 1.423e+02 3.211e+01 4.431 9.54e-06 ***
## hour8 4.603e+02 3.228e+01 14.260 < 2e-16 ***
## hour9 3.937e+01 3.319e+01 1.186 0.235628
## hour10 -1.783e+02 3.488e+01 -5.111 3.30e-07 ***
## hour11 -1.888e+02 3.606e+01 -5.235 1.70e-07 ***
## hour12 -1.321e+02 3.691e+01 -3.580 0.000346 ***
## hour13 -1.439e+02 3.776e+01 -3.813 0.000139 ***
## hour14 -1.487e+02 3.637e+01 -4.088 4.41e-05 ***
## hour15 -5.469e+01 3.586e+01 -1.525 0.127300
## hour16 1.149e+02 3.464e+01 3.318 0.000911 ***
## hour17 3.633e+02 3.347e+01 10.855 < 2e-16 ***
## hour18 8.332e+02 3.277e+01 25.427 < 2e-16 ***
## hour19 5.633e+02 3.269e+01 17.233 < 2e-16 ***
## hour20 4.896e+02 3.239e+01 15.115 < 2e-16 ***
## hour21 4.643e+02 3.225e+01 14.398 < 2e-16 ***
## hour22 3.768e+02 3.208e+01 11.746 < 2e-16 ***
## hour23 1.262e+02 3.186e+01 3.959 7.60e-05 ***
## temperature_c 3.806e+00 4.736e+00 0.804 0.421647
## humidity_percent -1.010e+01 1.347e+00 -7.499 7.32e-14 ***
## wind_speed_m_s 5.202e+00 5.424e+00 0.959 0.337584
## visibility_10m -3.460e-03 9.899e-03 -0.350 0.726679
## dew_point_temperature_c 2.007e+01 4.986e+00 4.025 5.76e-05 ***
## solar_radiation_mj_m2 6.831e+01 1.147e+01 5.955 2.74e-09 ***
## rainfall_mm -2.526e+02 1.956e+01 -12.917 < 2e-16 ***
## snowfall_cm 6.868e+01 1.945e+01 3.531 0.000418 ***
## seasonsSpring -1.824e+02 1.383e+01 -13.188 < 2e-16 ***
## seasonsSummer -1.750e+02 1.724e+01 -10.154 < 2e-16 ***
## seasonsWinter -4.089e+02 1.969e+01 -20.770 < 2e-16 ***
## holiday1 -1.356e+02 2.120e+01 -6.396 1.71e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 358.5 on 6261 degrees of freedom
## Multiple R-squared: 0.6835, Adjusted R-squared: 0.6818
## F-statistic: 386.4 on 35 and 6261 DF, p-value: < 2.2e-16Mô hình giải thích được hơn 68% lượng thay đổi trong các biến giải thích
anova(md.ln, test="Chisq")## Analysis of Variance Table
##
## Response: rented_bike_count
## Df Sum Sq Mean Sq F value Pr(>F)
## hour 23 852915014 37083261 288.5345 < 2.2e-16 ***
## temperature_c 1 699768135 699768135 5444.7008 < 2.2e-16 ***
## humidity_percent 1 62859972 62859972 489.0959 < 2.2e-16 ***
## wind_speed_m_s 1 83726 83726 0.6514 0.4196257
## visibility_10m 1 1030162 1030162 8.0154 0.0046529 **
## dew_point_temperature_c 1 1934492 1934492 15.0517 0.0001057 ***
## solar_radiation_mj_m2 1 6674927 6674927 51.9357 6.409e-13 ***
## rainfall_mm 1 22859938 22859938 177.8668 < 2.2e-16 ***
## snowfall_cm 1 3532 3532 0.0275 0.8683299
## seasons 3 84629981 28209994 219.4941 < 2.2e-16 ***
## holiday 1 5257906 5257906 40.9103 1.709e-10 ***
## Residuals 6261 804681192 128523
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Các biến hour, temperature_c, humidity_percent, dew_point_temperature_c, solar_radiation_mj_m2, rainfall_mm, seasons, và holiday có ý nghĩa thống kê cao và ảnh hưởng lớn đến số lượng xe đạp thuê.
- Các biến wind_speed_m_s và snowfall_cm không có ý nghĩa thống kê trong mô hình này.
y_pre_train <- predict(md.ln,trainset,type = "response")
y_pre_test <- predict(md.ln,testset ,type = "response")
rmse_train <- mean((trainset$rented_bike_count-y_pre_train) ^2) |> sqrt()
rmse_test <- mean((testset$rented_bike_count-y_pre_test) ^2) |> sqrt()
print(paste("RMSE của trainset:",rmse_train))## [1] "RMSE của trainset: 357.474509243905"print(paste("RMSE của testset:",rmse_test))## [1] "RMSE của testset: 369.427222345817"md.ln.cv <- errorest(rented_bike_count ~ .,
data = data.model,model = lm,estimator="cv",
est.para=control.errorest(k=5, predictions = TRUE))
md.ln.cv##
## Call:
## errorest.data.frame(formula = rented_bike_count ~ ., data = data.model,
## model = lm, estimator = "cv", est.para = control.errorest(k = 5,
## predictions = TRUE))
##
## 5-fold cross-validation estimator of root mean squared error
##
## Root mean squared error: 361.5928Dựa vào Rmse trên test, train, cũng như kết quả k-fold cross-validation, thì kết quả mô hình không bị overfitting nhưng rmse vẫn còn cao.
set.seed(12)
rmse_l <- rep(1,ncol(data.numeric))
for( i in 1:ncol(data.numeric)){
data.pca <- as.data.frame(PCA(data.numeric |> select(-rented_bike_count), ncp =i, graph = F)$ind$coord)
data.pca$rented_bike_count <- data.numeric$rented_bike_count
md.ln.cv <- errorest(rented_bike_count ~ .,
data = data.pca,model = lm,estimator="cv",
est.para=control.errorest(k=5, predictions = TRUE))
rmse_l[i] <- md.ln.cv$error
}which.min(rmse_l)## [1] 38rmse_l## [1] 572.9073 522.7233 522.5773 509.0930 498.7240 461.3164 447.7140 447.0399
## [9] 420.6869 418.3386 416.4121 416.5305 415.0307 413.1027 412.1735 410.1889
## [17] 407.2417 402.2952 401.9710 402.2343 400.7602 401.0530 400.9924 400.4795
## [25] 400.7529 397.9674 397.5725 397.8474 376.2425 376.5101 367.1026 363.5321
## [33] 361.7508 361.5474 361.6642 361.4610 361.4258 361.4140 361.5543Dù có PCA thì kết quả vẫn không hiệu quả hơn là mấy, vậy có thể kết luận dữ liệu không bị đa công tuyến quá nhiều.
Poisson regression
md.ps <- glm(rented_bike_count ~ .,
data = trainset,family = poisson)
summary(md.ps)##
## Call:
## glm(formula = rented_bike_count ~ ., family = poisson, data = trainset)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 7.301e+00 1.337e-02 546.144 < 2e-16 ***
## hour1 -2.439e-01 3.965e-03 -61.524 < 2e-16 ***
## hour2 -5.692e-01 4.364e-03 -130.405 < 2e-16 ***
## hour3 -9.130e-01 4.967e-03 -183.813 < 2e-16 ***
## hour4 -1.360e+00 5.805e-03 -234.262 < 2e-16 ***
## hour5 -1.283e+00 5.814e-03 -220.752 < 2e-16 ***
## hour6 -5.431e-01 4.453e-03 -121.955 < 2e-16 ***
## hour7 2.108e-01 3.615e-03 58.324 < 2e-16 ***
## hour8 5.965e-01 3.348e-03 178.167 < 2e-16 ***
## hour9 1.125e-01 3.728e-03 30.174 < 2e-16 ***
## hour10 -1.572e-01 4.063e-03 -38.682 < 2e-16 ***
## hour11 -1.179e-01 4.124e-03 -28.587 < 2e-16 ***
## hour12 -4.731e-03 4.107e-03 -1.152 0.24933
## hour13 -1.286e-02 4.156e-03 -3.095 0.00197 **
## hour14 -1.752e-02 4.029e-03 -4.348 1.37e-05 ***
## hour15 8.817e-02 3.894e-03 22.645 < 2e-16 ***
## hour16 2.524e-01 3.652e-03 69.119 < 2e-16 ***
## hour17 4.728e-01 3.424e-03 138.077 < 2e-16 ***
## hour18 8.230e-01 3.200e-03 257.208 < 2e-16 ***
## hour19 6.486e-01 3.257e-03 199.117 < 2e-16 ***
## hour20 5.929e-01 3.272e-03 181.198 < 2e-16 ***
## hour21 5.860e-01 3.283e-03 178.504 < 2e-16 ***
## hour22 4.884e-01 3.315e-03 147.314 < 2e-16 ***
## hour23 1.975e-01 3.523e-03 56.046 < 2e-16 ***
## temperature_c 1.589e-03 4.956e-04 3.206 0.00134 **
## humidity_percent -1.223e-02 1.496e-04 -81.762 < 2e-16 ***
## wind_speed_m_s -1.028e-02 5.772e-04 -17.805 < 2e-16 ***
## visibility_10m -2.278e-05 1.044e-06 -21.820 < 2e-16 ***
## dew_point_temperature_c 2.635e-02 5.240e-04 50.292 < 2e-16 ***
## solar_radiation_mj_m2 3.165e-02 1.171e-03 27.021 < 2e-16 ***
## rainfall_mm -1.246e+00 6.741e-03 -184.888 < 2e-16 ***
## snowfall_cm -7.930e-02 4.059e-03 -19.537 < 2e-16 ***
## seasonsSpring -2.068e-01 1.351e-03 -153.076 < 2e-16 ***
## seasonsSummer -2.033e-01 1.601e-03 -126.974 < 2e-16 ***
## seasonsWinter -1.019e+00 2.433e-03 -418.852 < 2e-16 ***
## holiday1 -1.812e-01 2.545e-03 -71.194 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 3247619 on 6296 degrees of freedom
## Residual deviance: 640109 on 6261 degrees of freedom
## AIC: 690550
##
## Number of Fisher Scoring iterations: 5Dựa vào bảng thống kê trên ta nhận thấy tham số phân tán là 1, nhưng thực chất với dữ liệu thực tế tình trạng này là thường xuyên gặp. Cụ thể trong trường hợp này là hơn 100 (\(\phi = \frac{\sum (\text{Pearson residuals})^2}{n - p}\)), do đấy mô hình đang bị overdispersion. Khi đó giả định về giả định rằng phương sai bằng trung bình không còn đúng.
anova(md.ps, test="Chisq")## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: rented_bike_count
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 6296 3247619
## hour 23 1248049 6273 1999569 < 2.2e-16 ***
## temperature_c 1 963614 6272 1035955 < 2.2e-16 ***
## humidity_percent 1 46470 6271 989485 < 2.2e-16 ***
## wind_speed_m_s 1 2585 6270 986900 < 2.2e-16 ***
## visibility_10m 1 234 6269 986666 < 2.2e-16 ***
## dew_point_temperature_c 1 4743 6268 981923 < 2.2e-16 ***
## solar_radiation_mj_m2 1 1657 6267 980266 < 2.2e-16 ***
## rainfall_mm 1 72578 6266 907688 < 2.2e-16 ***
## snowfall_cm 1 12412 6265 895276 < 2.2e-16 ***
## seasons 3 249812 6262 645464 < 2.2e-16 ***
## holiday 1 5355 6261 640109 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Các biến giải thích đều có p-value < 2.2e-16.
Biến giờ (hour), nhiệt độ (temperature_c), độ ẩm (humidity_percent), mùa (seasons) và lượng mưa (rainfall_mm) là những biến có ảnh hưởng lớn nhất.
Các biến khác như tốc độ gió, tầm nhìn, nhiệt độ điểm sương, bức xạ mặt trời, lượng tuyết rơi, và ngày lễ cũng có ảnh hưởng nhưng mức độ ít hơn so với các biến chính.
y_pre_train <- predict(md.ps,trainset,type = "response")
y_pre_test <- predict(md.ps,testset ,type = "response")
rmse_train <- mean((trainset$rented_bike_count-y_pre_train) ^2) |> sqrt()
rmse_test <- mean((testset$rented_bike_count-y_pre_test) ^2) |> sqrt()
print(paste("Rmse của trainset:",rmse_train))## [1] "Rmse của trainset: 312.076793691618"print(paste("Rmse của testset:",rmse_test))## [1] "Rmse của testset: 320.641281267276"myGLM = function(formula, data) {
glm(formula, data, family = poisson(link = log))
}
myPredictGLM = function(object, newdata){
predict(object, newdata , type="response")
}
md.ps.cv <- errorest(rented_bike_count ~ .,
data = data.model, predict = myPredictGLM,model = myGLM,estimator="cv",
est.para=control.errorest(k=5, predictions = TRUE))
md.ps.cv##
## Call:
## errorest.data.frame(formula = rented_bike_count ~ ., data = data.model,
## model = myGLM, predict = myPredictGLM, estimator = "cv",
## est.para = control.errorest(k = 5, predictions = TRUE))
##
## 5-fold cross-validation estimator of root mean squared error
##
## Root mean squared error: 315.2939Kết quả rmse tốt hơn so với mô hình hồi quy tuyến tính.
md.ps.null <- glm(rented_bike_count ~ 1, family = "poisson", data = trainset)
md.ps.R2 <- 1 - as.numeric(logLik(md.ps)) / as.numeric(logLik(md.ps.null))
md.ps.R2## [1] 0.7906365Với \(Pseudo-R^2\) gần 0.8 cho thấy mức giải thích được cải thiện là gần 80% với vo mô hình chỉ có bias.
Quasipoisson
md.qps <- glm(rented_bike_count ~ .,
data = trainset,family = quasipoisson)
summary(md.qps)##
## Call:
## glm(formula = rented_bike_count ~ ., family = quasipoisson, data = trainset)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.301e+00 1.348e-01 54.148 < 2e-16 ***
## hour1 -2.439e-01 3.999e-02 -6.100 1.13e-09 ***
## hour2 -5.692e-01 4.402e-02 -12.929 < 2e-16 ***
## hour3 -9.130e-01 5.010e-02 -18.224 < 2e-16 ***
## hour4 -1.360e+00 5.856e-02 -23.226 < 2e-16 ***
## hour5 -1.283e+00 5.864e-02 -21.887 < 2e-16 ***
## hour6 -5.431e-01 4.492e-02 -12.091 < 2e-16 ***
## hour7 2.108e-01 3.646e-02 5.783 7.71e-09 ***
## hour8 5.965e-01 3.377e-02 17.664 < 2e-16 ***
## hour9 1.125e-01 3.760e-02 2.992 0.002786 **
## hour10 -1.572e-01 4.098e-02 -3.835 0.000127 ***
## hour11 -1.179e-01 4.159e-02 -2.834 0.004608 **
## hour12 -4.731e-03 4.142e-02 -0.114 0.909073
## hour13 -1.286e-02 4.192e-02 -0.307 0.758998
## hour14 -1.752e-02 4.064e-02 -0.431 0.666429
## hour15 8.817e-02 3.927e-02 2.245 0.024795 *
## hour16 2.524e-01 3.683e-02 6.853 7.93e-12 ***
## hour17 4.728e-01 3.453e-02 13.690 < 2e-16 ***
## hour18 8.230e-01 3.227e-02 25.501 < 2e-16 ***
## hour19 6.486e-01 3.285e-02 19.741 < 2e-16 ***
## hour20 5.929e-01 3.300e-02 17.965 < 2e-16 ***
## hour21 5.860e-01 3.311e-02 17.698 < 2e-16 ***
## hour22 4.884e-01 3.344e-02 14.605 < 2e-16 ***
## hour23 1.975e-01 3.553e-02 5.557 2.86e-08 ***
## temperature_c 1.589e-03 4.999e-03 0.318 0.750575
## humidity_percent -1.223e-02 1.509e-03 -8.106 6.23e-16 ***
## wind_speed_m_s -1.028e-02 5.822e-03 -1.765 0.077573 .
## visibility_10m -2.278e-05 1.053e-05 -2.163 0.030552 *
## dew_point_temperature_c 2.635e-02 5.285e-03 4.986 6.32e-07 ***
## solar_radiation_mj_m2 3.165e-02 1.181e-02 2.679 0.007403 **
## rainfall_mm -1.246e+00 6.800e-02 -18.331 < 2e-16 ***
## snowfall_cm -7.930e-02 4.094e-02 -1.937 0.052786 .
## seasonsSpring -2.068e-01 1.362e-02 -15.177 < 2e-16 ***
## seasonsSummer -2.033e-01 1.615e-02 -12.589 < 2e-16 ***
## seasonsWinter -1.019e+00 2.454e-02 -41.527 < 2e-16 ***
## holiday1 -1.812e-01 2.567e-02 -7.059 1.86e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 101.7313)
##
## Null deviance: 3247619 on 6296 degrees of freedom
## Residual deviance: 640109 on 6261 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 5anova(md.qps, test="Chisq")## Analysis of Deviance Table
##
## Model: quasipoisson, link: log
##
## Response: rented_bike_count
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 6296 3247619
## hour 23 1248049 6273 1999569 < 2.2e-16 ***
## temperature_c 1 963614 6272 1035955 < 2.2e-16 ***
## humidity_percent 1 46470 6271 989485 < 2.2e-16 ***
## wind_speed_m_s 1 2585 6270 986900 4.640e-07 ***
## visibility_10m 1 234 6269 986666 0.1293
## dew_point_temperature_c 1 4743 6268 981923 8.591e-12 ***
## solar_radiation_mj_m2 1 1657 6267 980266 5.449e-05 ***
## rainfall_mm 1 72578 6266 907688 < 2.2e-16 ***
## snowfall_cm 1 12412 6265 895276 < 2.2e-16 ***
## seasons 3 249812 6262 645464 < 2.2e-16 ***
## holiday 1 5355 6261 640109 4.009e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Các biến hour, temperature_c, humidity_percent, wind_speed_m_s, dew_point_temperature_c, solar_radiation_mj_m2, rainfall_mm, snowfall_cm, seasons, và holiday đều có mức ý nghĩa thống kê cao (Pr(>Chi) < 0.05).
- Biến visibility_10m không có mức ý nghĩa thống kê (Pr(>Chi) > 0.05), cho thấy nó không có ảnh hưởng đáng kể đến số lượng xe đạp thuê trong mô hình này.
- Dữ liệu còn khá phân tán có thể chia ra thành các mùa, để các dữ liệu ít bị phân tác hơn.
y_pre_train <- predict(md.qps,trainset,type = "response")
y_pre_test <- predict(md.qps,testset ,type = "response")
rmse_train <- mean((trainset$rented_bike_count-y_pre_train) ^2) |> sqrt()
rmse_test <- mean((testset$rented_bike_count-y_pre_test) ^2) |> sqrt()
print(paste("Rmse của trainset:",rmse_train))## [1] "Rmse của trainset: 312.076793691618"print(paste("Rmse của testset:",rmse_test))## [1] "Rmse của testset: 320.641281267276"Đánh giá
autoplot(md.ps, which = 1, ncol = 1, label.size = 3,
colour = "seasons",alpha = 0.6) + theme_bw()
Các điểm thặng dư đối xứng qua trục y=0, không xuất hiện đường cong bất thường nào. Điều này cho thấy tính tuyến tính của mô hình được bảo đảm. Tuy nhiên giả định đồng nhất phương sai là không đảm bảo, vì thặng dư co cụm về bên phải tăng đều, và phân tán mạnh hơn. Điều này một phần vì dữ liệu phân tán khá lớn, tuy nhiên các mô hình cũng không cần giả định này phải thoả mãn
autoplot(md.ps, which = 2, ncol = 1, label.size = 3,
colour = "seasons",alpha = 0.6) + theme_bw()
Thặng dư không tuân theo phân phối chuẩn, điều này là hiển nhiên vì đây là hồi quy poisson, nhưng số lượng mẫu khá lớn nên thăng dư cũng xấp xỉ gần với phân phối chuẩn.
autoplot(md.ps, which = 3, ncol = 1, label.size = 3,
colour = "seasons",alpha = 0.6) + theme_bw()
Như đã đề cập ở dữ liệu phân tán khá cao, giả định về đồng nhất phương sai thặng dư không đảm bảo. Ở đây ta thấy rõ ràng ngoài xu hướng càng tăng giá trị fit, càng tăng dần độ phân tán, độ phân tán còn có sự khác nhau ở các mùa.
autoplot(md.ps, which = 4:6, ncol = 1, label.size = 3,
colour = "seasons",alpha = 0.6) + theme_bw()
Cook’s Distance: dùng để đo lường mức độ ảnh hưởng của từng điểm dữ liệu đối với mô hình hồi quy. Điểm dữ liệu có giá trị Cook’s Distance lớn có thể có ảnh hưởng lớn đến các tham số ước lượng của mô hình. +) Điểm 5212, 4494, 6029: Các điểm này có giá trị Cook’s Distance cao hơn hẳn so với các điểm khác, cho thấy chúng có ảnh hưởng lớn đến mô hình. Những điểm này có thể là outliers hoặc có leverage cao. +) Phần lớn các điểm có giá trị Cook’s Distance nhỏ, cho thấy chúng không ảnh hưởng nhiều đến mô hình.
Residuals vs Leverage: kiểm tra mối quan hệ giữa residuals (sai số dự đoán) và leverage (mức độ ảnh hưởng của một quan sát trong việc xác định ước lượng của mô hình). +) Residuals lớn (Điểm 5212, 6029, 4494): Các điểm này có residuals lớn, có thể là dấu hiệu của các outliers, cho thấy mô hình dự đoán kém đối với những điểm này. +) Phân bố residuals: Phân bố residuals không đối xứng hoàn toàn xung quanh trục y = 0. Điều này có thể chỉ ra rằng mô hình có vấn đề về phân phối sai số. +) Leverage thấp: Hầu hết các điểm có leverage thấp (gần 0), cho thấy chúng không ảnh hưởng mạnh đến mô hình.
Cook’s Distance vs Leverage: kiểm tra mối quan hệ giữa Cook’s Distance và leverage. Lệu những điểm có leverage cao (ảnh hưởng lớn) cũng có Cook’s Distance cao (ảnh hưởng đến mô hình) hay không. +) Điểm 5212, 6029, 4494: Các điểm này có cả Cook’s Distance và leverage cao, cho thấy chúng có ảnh hưởng lớn đến mô hình và cần được xem xét kỹ lưỡng. Nếu các điểm này là các outliers, có thể cân nhắc loại bỏ chúng hoặc kiểm tra kỹ hơn để xem chúng có hợp lý hay không. +) Phần lớn các điểm: Hầu hết các điểm còn lại có Cook’s Distance và leverage thấp, cho thấy chúng không ảnh hưởng nhiều đến mô hình.
data.md.vs <- data.frame(trainset |> select(temperature_c,humidity_percent,wind_speed_m_s,visibility_10m,rainfall_mm,snowfall_cm,seasons) , parital_residual = residuals(md.ps) )
data.md.vs.m <- gather(data.md.vs,var,val,-c(seasons,parital_residual) )
ggplot(data.md.vs.m,aes(x = val, y = parital_residual, color = seasons))+facet_wrap(~var,scales = 'free_x') + geom_point()
Humidity Percent +) Các giá trị partial residuals phân bố khá đồng đều quanh trục y = 0. +) Không có sự khác biệt rõ ràng về phân bố partial residuals theo mùa. +) Độ ẩm có vẻ không có mối quan hệ phi tuyến với biến đáp ứng. Có thể biến này được mô hình hóa tốt bởi Poisson Regression.
Rainfall mm +) Các giá trị partial residuals khá phân tán, đặc biệt ở mức lượng mưa thấp (dưới 1 mm). +) Phân bố partial residuals không thay đổi nhiều theo mùa. +) Mối quan hệ không rõ ràng, có thể do một số outliers. Cần kiểm tra thêm hoặc có thể cân nhắc biến đổi biến này
Snowfall cm +) Phân bố partial residuals: Các giá trị partial residuals chủ yếu tập trung ở mức tuyết rơi bằng 0, với một số outliers ở các giá trị tuyết rơi cao hơn. +) Phần lớn các điểm tuyết rơi lớn xuất hiện vào mùa đông. +) Mối quan hệ phi tuyến rõ ràng hơn, đặc biệt vào mùa đông. Cần kiểm tra thêm và có thể cần một số biến đổi hoặc thêm biến tương tác để cải thiện mô hình.
Temperature c +) Các giá trị partial residuals phân bố khá đồng đều quanh trục y = 0. +) Có sự phân bố rõ ràng theo mùa: nhiệt độ cao hơn vào mùa hè và thấp hơn vào mùa đông. +) Mối quan hệ giữa nhiệt độ và partial residuals có thể được biểu diễn bằng một mô hình tuyến tính.
Visibility 10m +) Các giá trị partial residuals phân bố đồng đều, với một số outliers. +) Không có sự khác biệt rõ ràng về phân bố partial residuals theo mùa. +) Mối quan hệ không rõ ràng. Cần kiểm tra thêm để xem liệu biến này có cần biến đổi hoặc loại bỏ không.
Wind Speed m/s +) Các giá trị partial residuals phân bố đồng đều quanh trục y = 0. +) Không có sự khác biệt rõ ràng về phân bố partial residuals theo mùa. +) Tốc độ gió có vẻ được mô hình hóa tốt bởi Poisson Regression. Mối quan hệ tuyến tính có thể hợp lý.
Khoảng tin cậy bootstraping:
fun.boot<- function(data, ind, formula, ...){
data_new <- data[ind,]
out_md <- glm(formula = formula, data = data_new,family = poisson,...)
return(out_md$coefficients)
}
md.ps.boot <- boot(data = data.model, statistic = fun.boot, R = 500,
formula = rented_bike_count~.)
summary(md.ps.boot)##
## Number of bootstrap replications R = 500
## original bootBias bootSE bootMed
## 1 7.2013e+00 -3.7976e-03 1.3985e-01 7.1927e+00
## 2 -2.4136e-01 2.8829e-03 2.5424e-02 -2.3931e-01
## 3 -5.7418e-01 1.4716e-03 2.7498e-02 -5.7214e-01
## 4 -9.2478e-01 2.4856e-04 3.0532e-02 -9.2251e-01
## 5 -1.3666e+00 1.9385e-03 2.7247e-02 -1.3637e+00
## 6 -1.2869e+00 3.2671e-04 2.8212e-02 -1.2847e+00
## 7 -5.3392e-01 2.0532e-03 3.2465e-02 -5.3123e-01
## 8 1.9199e-01 2.9494e-03 3.4654e-02 1.9531e-01
## 9 6.3545e-01 2.3324e-03 3.5198e-02 6.3773e-01
## 10 1.1536e-01 1.5840e-03 2.7225e-02 1.1463e-01
## 11 -1.7923e-01 2.5840e-03 2.9299e-02 -1.7710e-01
## 12 -1.1787e-01 2.3996e-03 3.1779e-02 -1.1614e-01
## 13 -1.2281e-02 1.8045e-03 3.4071e-02 -1.0055e-02
## 14 -2.6172e-02 1.3404e-03 3.5162e-02 -2.4951e-02
## 15 -2.1157e-02 2.7690e-03 3.3273e-02 -1.8601e-02
## 16 7.7703e-02 2.9788e-03 3.4197e-02 8.1346e-02
## 17 2.3722e-01 6.7044e-04 3.2202e-02 2.3882e-01
## 18 4.8965e-01 3.8637e-04 2.6837e-02 4.8986e-01
## 19 8.2282e-01 2.6467e-03 2.5710e-02 8.2455e-01
## 20 6.5895e-01 1.5563e-03 2.7639e-02 6.6107e-01
## 21 5.9000e-01 2.1887e-03 2.5735e-02 5.9137e-01
## 22 5.8623e-01 2.3546e-03 2.5011e-02 5.8968e-01
## 23 4.8602e-01 2.3758e-03 2.6354e-02 4.8956e-01
## 24 1.9370e-01 1.7097e-03 2.6574e-02 1.9476e-01
## 25 4.1078e-03 9.1163e-05 5.0174e-03 4.2269e-03
## 26 -1.1191e-02 3.0218e-05 1.6019e-03 -1.1167e-02
## 27 -8.9989e-03 -1.2178e-04 5.5612e-03 -9.0632e-03
## 28 -1.3289e-05 -4.8165e-08 1.0229e-05 -1.3441e-05
## 29 2.3357e-02 -8.0026e-05 5.2639e-03 2.3203e-02
## 30 3.3739e-02 -3.1836e-04 1.1250e-02 3.3410e-02
## 31 -1.3279e+00 -7.2986e-03 1.2356e-01 -1.3333e+00
## 32 -8.1502e-02 -4.2052e-04 2.0383e-02 -8.0285e-02
## 33 -1.9582e-01 -1.1999e-03 1.3003e-02 -1.9600e-01
## 34 -1.9601e-01 -1.7357e-04 1.5688e-02 -1.9461e-01
## 35 -1.0123e+00 4.2797e-04 2.2067e-02 -1.0117e+00
## 36 -1.7515e-01 -1.9252e-04 3.0740e-02 -1.7603e-01boot.ci(md.ps.boot, index = 1, type = "perc", conf = 0.95)## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 500 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = md.ps.boot, conf = 0.95, type = "perc", index = 1)
##
## Intervals :
## Level Percentile
## 95% ( 6.939, 7.495 )
## Calculations and Intervals on Original Scaleset.seed(14)
sample <- data.model |> select(-rented_bike_count) |> sample_n(1) |> mutate(across(where(is.numeric), ~ . + .*rnorm(1, mean = 0, sd = 0.01)))
print(t(sample))## [,1]
## hour "21"
## temperature_c "-8.831034"
## humidity_percent "43.0568"
## wind_speed_m_s "1.400403"
## visibility_10m "1995.394"
## dew_point_temperature_c "-19.04248"
## solar_radiation_mj_m2 "0"
## rainfall_mm "0"
## snowfall_cm "0"
## seasons "Winter"
## holiday "0"sample <- model.matrix(~ ., data = sample)
y.pre.boot <- apply(md.ps.boot$t, 1, function(x){
exp(sample %*% x)
})
quantile(y.pre.boot, probs = c(0.025, 0.975))## 2.5% 97.5%
## 307.7393 335.3752Cải thiện model
Mô hình cho từng mùa
Winter
data.model.wt <- data.model[ data.model$ seasons =="Winter",- which(colnames(data.model) %in% c('seasons'))]
rs <- slit_test_train(data.model.wt,0.8)
trainset <- as.data.frame(rs[1])
testset <- as.data.frame(rs[2])
md.qps.wt <- glm(rented_bike_count ~ .,family = quasipoisson('log'),data=trainset )
summary(md.qps.wt)##
## Call:
## glm(formula = rented_bike_count ~ ., family = quasipoisson("log"),
## data = trainset)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.102e+00 2.293e-01 26.610 < 2e-16 ***
## hour1 -7.128e-03 6.951e-02 -0.103 0.918331
## hour2 -3.884e-01 7.420e-02 -5.234 1.87e-07 ***
## hour3 -7.416e-01 8.407e-02 -8.822 < 2e-16 ***
## hour4 -1.167e+00 9.796e-02 -11.909 < 2e-16 ***
## hour5 -1.132e+00 9.463e-02 -11.958 < 2e-16 ***
## hour6 -4.981e-01 7.681e-02 -6.485 1.18e-10 ***
## hour7 3.253e-01 6.351e-02 5.121 3.41e-07 ***
## hour8 9.490e-01 5.984e-02 15.858 < 2e-16 ***
## hour9 5.194e-01 6.270e-02 8.284 2.50e-16 ***
## hour10 1.298e-01 6.823e-02 1.902 0.057336 .
## hour11 2.351e-01 6.917e-02 3.398 0.000694 ***
## hour12 3.132e-01 7.061e-02 4.435 9.82e-06 ***
## hour13 3.348e-01 7.299e-02 4.587 4.85e-06 ***
## hour14 3.292e-01 6.978e-02 4.718 2.59e-06 ***
## hour15 3.527e-01 6.733e-02 5.238 1.83e-07 ***
## hour16 3.965e-01 6.548e-02 6.055 1.74e-09 ***
## hour17 5.296e-01 6.146e-02 8.618 < 2e-16 ***
## hour18 8.481e-01 5.839e-02 14.525 < 2e-16 ***
## hour19 4.945e-01 6.119e-02 8.081 1.26e-15 ***
## hour20 2.730e-01 6.418e-02 4.253 2.23e-05 ***
## hour21 3.054e-01 6.262e-02 4.878 1.18e-06 ***
## hour22 2.630e-01 6.477e-02 4.060 5.13e-05 ***
## hour23 3.175e-02 6.679e-02 0.475 0.634626
## temperature_c 1.591e-02 8.364e-03 1.903 0.057282 .
## humidity_percent -1.102e-02 2.589e-03 -4.256 2.20e-05 ***
## wind_speed_m_s -3.345e-02 8.443e-03 -3.962 7.77e-05 ***
## visibility_10m 1.625e-04 2.537e-05 6.405 1.98e-10 ***
## dew_point_temperature_c 3.819e-02 9.051e-03 4.219 2.59e-05 ***
## solar_radiation_mj_m2 -4.150e-03 3.073e-02 -0.135 0.892598
## rainfall_mm -3.069e-01 1.014e-01 -3.025 0.002527 **
## snowfall_cm -7.377e-02 2.196e-02 -3.360 0.000799 ***
## holiday1 -4.533e-01 3.432e-02 -13.210 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 24.48278)
##
## Null deviance: 147711 on 1627 degrees of freedom
## Residual deviance: 41915 on 1595 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 5y_pre_train <- predict(md.qps.wt,trainset,type = "response")
y_pre_test <- predict(md.qps.wt,testset ,type = "response")
rmse_train <- mean((trainset$rented_bike_count-y_pre_train) ^2) |> sqrt()
rmse_test <- mean((testset$rented_bike_count-y_pre_test) ^2) |> sqrt()
print(paste("Rmse của trainset:",rmse_train))## [1] "Rmse của trainset: 83.3299356536891"print(paste("Rmse của testset:",rmse_test))## [1] "Rmse của testset: 95.7229074446783"myGLM = function(formula, data) {
glm(formula, data, family = quasipoisson(link = log))
}
myPredictGLM = function(object, newdata){
predict(object, newdata , type="response")
}
md.qps.cv.wt <- errorest(rented_bike_count ~ .,
data = data.model.wt, predict = myPredictGLM,model = myGLM,estimator="cv",
est.para=control.errorest(k=5, predictions = TRUE))
md.qps.cv.wt##
## Call:
## errorest.data.frame(formula = rented_bike_count ~ ., data = data.model.wt,
## model = myGLM, predict = myPredictGLM, estimator = "cv",
## est.para = control.errorest(k = 5, predictions = TRUE))
##
## 5-fold cross-validation estimator of root mean squared error
##
## Root mean squared error: 87.8302Summer
data.model.sm <- data.model[ data.model$ seasons =="Summer",- which(colnames(data.model) %in% c('seasons'))]
rs <- slit_test_train(data.model.sm,0.8)
trainset <- as.data.frame(rs[1])
testset <- as.data.frame(rs[2])
md.qps.sm <- glm(rented_bike_count ~ .,family = quasipoisson('log'),data=trainset )
summary(md.qps.sm)##
## Call:
## glm(formula = rented_bike_count ~ ., family = quasipoisson("log"),
## data = trainset)
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.574e+00 2.559e-01 37.416 < 2e-16 ***
## hour1 -3.120e-01 4.996e-02 -6.246 5.41e-10 ***
## hour2 -6.352e-01 5.509e-02 -11.530 < 2e-16 ***
## hour3 -1.063e+00 6.423e-02 -16.548 < 2e-16 ***
## hour4 -1.468e+00 7.467e-02 -19.660 < 2e-16 ***
## hour5 -1.357e+00 7.392e-02 -18.361 < 2e-16 ***
## hour6 -7.258e-01 5.821e-02 -12.467 < 2e-16 ***
## hour7 -3.902e-02 4.723e-02 -0.826 0.4088
## hour8 2.510e-01 4.454e-02 5.636 2.06e-08 ***
## hour9 -1.987e-01 5.042e-02 -3.941 8.48e-05 ***
## hour10 -4.527e-01 5.529e-02 -8.187 5.43e-16 ***
## hour11 -3.784e-01 5.618e-02 -6.735 2.28e-11 ***
## hour12 -2.834e-01 5.707e-02 -4.966 7.59e-07 ***
## hour13 -2.695e-01 5.909e-02 -4.561 5.49e-06 ***
## hour14 -2.321e-01 5.683e-02 -4.085 4.64e-05 ***
## hour15 -8.945e-02 5.381e-02 -1.662 0.0967 .
## hour16 9.731e-02 5.075e-02 1.917 0.0554 .
## hour17 4.401e-01 4.573e-02 9.623 < 2e-16 ***
## hour18 8.409e-01 4.106e-02 20.479 < 2e-16 ***
## hour19 7.612e-01 4.097e-02 18.580 < 2e-16 ***
## hour20 7.137e-01 4.045e-02 17.647 < 2e-16 ***
## hour21 6.627e-01 4.043e-02 16.391 < 2e-16 ***
## hour22 5.279e-01 4.053e-02 13.026 < 2e-16 ***
## hour23 2.318e-01 4.323e-02 5.361 9.47e-08 ***
## temperature_c -8.505e-02 9.592e-03 -8.867 < 2e-16 ***
## humidity_percent -2.089e-02 2.856e-03 -7.314 4.09e-13 ***
## wind_speed_m_s 5.763e-03 8.675e-03 0.664 0.5066
## visibility_10m -8.665e-05 1.548e-05 -5.598 2.54e-08 ***
## dew_point_temperature_c 5.590e-02 1.015e-02 5.509 4.20e-08 ***
## solar_radiation_mj_m2 9.602e-02 1.674e-02 5.735 1.16e-08 ***
## rainfall_mm -1.075e+00 8.136e-02 -13.217 < 2e-16 ***
## snowfall_cm NA NA NA NA
## holiday1 -3.540e-02 4.355e-02 -0.813 0.4165
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 68.63507)
##
## Null deviance: 653959 on 1617 degrees of freedom
## Residual deviance: 116143 on 1586 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 5y_pre_train <- predict(md.qps.sm,trainset,type = "response")
y_pre_test <- predict(md.qps.sm,testset ,type = "response")
rmse_train <- mean((trainset$rented_bike_count-y_pre_train) ^2) |> sqrt()
rmse_test <- mean((testset$rented_bike_count-y_pre_test) ^2) |> sqrt()
print(paste("Rmse của trainset:",rmse_train))## [1] "Rmse của trainset: 287.701859290571"print(paste("Rmse của testset:",rmse_test))## [1] "Rmse của testset: 320.851827259716"myGLM = function(formula, data) {
glm(formula, data, family = quasipoisson(link = log))
}
myPredictGLM = function(object, newdata){
predict(object, newdata , type="response")
}
md.qps.cv.sm <- errorest(rented_bike_count ~ .,
data = data.model.sm, predict = myPredictGLM,model = myGLM,estimator="cv",
est.para=control.errorest(k=5, predictions = TRUE))
md.qps.cv.sm##
## Call:
## errorest.data.frame(formula = rented_bike_count ~ ., data = data.model.sm,
## model = myGLM, predict = myPredictGLM, estimator = "cv",
## est.para = control.errorest(k = 5, predictions = TRUE))
##
## 5-fold cross-validation estimator of root mean squared error
##
## Root mean squared error: 298.6332Sau khi phân chia ra các mùa cụ thể ta thấy dữ liệu trong từng nhóm đồng chất hơn nên độ phân tán mô hình thấp hơn và cũng như rmse cũng thấp mô hình tổng.
Splines
Cải thiện cho biến humidity_percent và temperature_c là 2 biến tác động lớn trong mô hình hồi quy poisson
rs <- slit_test_train(data.model,0.8)
trainset <- as.data.frame(rs[1])
testset <- as.data.frame(rs[2])
knots_h <- quantile(trainset$humidity_percent, probs = c(0.25, 0.5, 0.75))
knots_t <- quantile(trainset$temperature_c, probs = c(0.25, 0.5, 0.75))
md.ps.sl <- glm(rented_bike_count ~ hour + bs(temperature_c,knots_t,3) + bs(humidity_percent,knots_h,3) + wind_speed_m_s + dew_point_temperature_c + solar_radiation_mj_m2 + rainfall_mm + snowfall_cm + holiday + seasons,
data = trainset,family = quasipoisson)
y_pre_train <- predict(md.ps.sl,trainset,type = "response")
y_pre_test <- predict(md.ps.sl,testset ,type = "response")
rmse_train <- mean((trainset$rented_bike_count-y_pre_train) ^2) |> sqrt()
rmse_test <- mean((testset$rented_bike_count-y_pre_test) ^2) |> sqrt()
print(paste("Rmse của trainset:",rmse_train))## [1] "Rmse của trainset: 253.687888119918"print(paste("Rmse của testset:",rmse_test))## [1] "Rmse của testset: 259.658732815563"print(paste("Độ phân tán:",sum(residuals(md.ps.sl,"pearson" ) ^2 ) /md.ps.sl$df.residual ))## [1] "Độ phân tán: 73.9049528539409"Mô hình có sự tương tác
rs <- slit_test_train(data.model,0.8)
trainset <- as.data.frame(rs[1])
testset <- as.data.frame(rs[2])
md.ps.int <- glm(rented_bike_count ~ (hour + temperature_c + humidity_percent + wind_speed_m_s + dew_point_temperature_c + solar_radiation_mj_m2 + rainfall_mm + snowfall_cm + holiday) * seasons,
data = trainset,family = quasipoisson)
y_pre_train <- predict(md.ps.int,trainset,type = "response")
y_pre_test <- predict(md.ps.int,testset ,type = "response")
rmse_train <- mean((trainset$rented_bike_count-y_pre_train) ^2) |> sqrt()
rmse_test <- mean((testset$rented_bike_count-y_pre_test) ^2) |> sqrt()
print(paste("Rmse của trainset:",rmse_train))## [1] "Rmse của trainset: 247.393906091801"print(paste("Rmse của testset:",rmse_test))## [1] "Rmse của testset: 253.290239415538"print(paste("Độ phân tán:",sum(residuals(md.ps.int,"pearson" ) ^2 ) /md.ps.int$df.residual ))## [1] "Độ phân tán: 70.4196996642662"anova(md.ps.int, test="Chisq")## Analysis of Deviance Table
##
## Model: quasipoisson, link: log
##
## Response: rented_bike_count
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 6296 3247619
## hour 23 1248049 6273 1999569 < 2.2e-16 ***
## temperature_c 1 963614 6272 1035955 < 2.2e-16 ***
## humidity_percent 1 46470 6271 989485 < 2.2e-16 ***
## wind_speed_m_s 1 2585 6270 986900 1.380e-09 ***
## dew_point_temperature_c 1 4947 6269 981953 < 2.2e-16 ***
## solar_radiation_mj_m2 1 1586 6268 980367 2.076e-06 ***
## rainfall_mm 1 72584 6267 907783 < 2.2e-16 ***
## snowfall_cm 1 12329 6266 895454 < 2.2e-16 ***
## holiday 1 9547 6265 885907 < 2.2e-16 ***
## seasons 3 245323 6262 640584 < 2.2e-16 ***
## hour:seasons 69 75782 6193 564802 < 2.2e-16 ***
## temperature_c:seasons 3 122509 6190 442293 < 2.2e-16 ***
## humidity_percent:seasons 3 3311 6187 438982 3.468e-10 ***
## wind_speed_m_s:seasons 3 1568 6184 437414 5.748e-05 ***
## dew_point_temperature_c:seasons 3 1461 6181 435953 0.0001194 ***
## solar_radiation_mj_m2:seasons 3 480 6178 435473 0.0778307 .
## rainfall_mm:seasons 3 470 6175 435003 0.0828562 .
## snowfall_cm:seasons 1 936 6174 434066 0.0002666 ***
## holiday:seasons 3 3235 6171 430831 5.878e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Mô hình khi có sự tương tác giữa seasons với các yếu tố còn lại thực sự hiểu quả hơn, tuy nhiên độ phức tạp cũng tăng.
Co hệ số
X <- data.numeric |> select(-rented_bike_count) |> scale() |> as.matrix()
y <- data.numeric |> select(rented_bike_count) |> unlist()
out.cv.lasso <- cv.glmnet(x = X, y = y, alpha = 1,
type.measure = "mse", nfolds = 5,
family = poisson)
print(out.cv.lasso)##
## Call: cv.glmnet(x = X, y = y, type.measure = "mse", nfolds = 5, alpha = 1, family = poisson)
##
## Measure: Mean-Squared Error
##
## Lambda Index Measure SE Nonzero
## min 0.421 74 99637 2447 35
## 1se 6.257 45 101888 2333 29beta.lambda.lasso <- out.cv.lasso$lambda.min
md.lasso <- glmnet(x = X, y = y, alpha = 1, family = quasipoisson)
predict(md.lasso, s = beta.lambda.lasso, type = "coefficients")## 39 x 1 sparse Matrix of class "dgCMatrix"
## s1
## (Intercept) 6.261395e+00
## hour.0 -1.526357e-02
## hour.1 -6.417294e-02
## hour.2 -1.293567e-01
## hour.3 -2.001664e-01
## hour.4 -2.892629e-01
## hour.5 -2.706017e-01
## hour.6 -1.218734e-01
## hour.7 2.201748e-02
## hour.8 1.125281e-01
## hour.9 7.805370e-03
## hour.10 -4.903058e-02
## hour.11 -3.680459e-02
## hour.12 -1.586947e-02
## hour.13 -1.882491e-02
## hour.14 -1.817354e-02
## hour.15 .
## hour.16 3.154337e-02
## hour.17 8.232566e-02
## hour.18 1.489885e-01
## hour.19 1.145066e-01
## hour.20 1.006405e-01
## hour.21 1.008192e-01
## hour.22 8.100793e-02
## hour.23 2.202673e-02
## temperature_c 1.116877e-01
## humidity_percent -1.881183e-01
## wind_speed_m_s -7.145285e-03
## visibility_10m -5.719105e-03
## dew_point_temperature_c 2.401607e-01
## solar_radiation_mj_m2 2.427883e-02
## rainfall_mm -3.306467e-01
## snowfall_cm -1.845556e-02
## seasons.Autumn 8.104036e-02
## seasons.Spring .
## seasons.Summer .
## seasons.Winter -3.564635e-01
## holiday.0 3.714379e-02
## holiday.1 -6.736188e-10Với phương pháp co hệ số ta co thể loại bỏ một số biến không tác động đến mô hình.
Tổng kết:
Sau thực hiện mô hình ta rút ra được những insight, đó là sự biến đổi số lượng xe thuê giá phụ thuộc vào mùa, thời gian,và sự tác động các biến giải thích đến số lượng xe thuê. Từ đấy, phục vụ cho các mục đích kinh doanh, marketing,…